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Abstract 

High multiplicity events in p+p collisions are studied using the theory of the Color Glass Condensate. We 
show that intrinsic fluctuations of the proton saturation momentum scale are needed in addition to the 
sub-nucleonic color charge fluctuations to explain the very high multiplicity tail of distributions in p+p 
collisions. The origin of such intrinsic fluctuations are presumably non-perturbative in nature. Classical 
Yang Mills simulations using the IP-Glasma model are performed to make quantitative estimations. We 
find that fluctuations as large as Oil) of the average values of the saturation momentum scale can lead 
to rare high multiplicity events seen in p+p data at RHIC and LHC energies. Using the available data on 
multiplicity distributions we try to constrain the distribution of the proton saturation momentum scale 
and make predictions for the multiplicity distribution in 13 TeV p+p collisions. 


1. Introduction 

High multiplicity events in the collisions of small systems like p+p have recently become an important 
topic of interest. Such processes can provide important insights about the dynamics of multi-particle 
production arising from the saturated gluon fields inside the protons. However, the origin of such high 
multiplicity events from the first principles QCD approach is not yet fully understood. We demonstrate 
in this paper that the description of multiplicity fluctuations requires a framework of multi-particle 
production that includes several non-perturbative scales in the problem. 

Recent efforts to understand multpilicity fluctuations have focused on applying the theory of the 
Color Glass Condensate (CGC) [TJ El [3l IH [5]. This theory successfully describes several features of the 
global data in hadronic, light-heavy, and heavy ion collisions. In the CGC framework, the emergence of 
a dynamical scale Qs makes the computation of multi-particle production feasible at very high energies. 
This is because the large values of Qs make the effective coupling small, thereby making weak coupling 
methods applicable. Qs being the only scale in the problem also naturally sets the scale of gluon field 
fluctuations that eventually give rise to the fluctuations in multiplicity. It turns out that using a single 
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framework one can qualitatively describe both the bulk multiplicity and its fluctuations up to the level of 
uncertainties introduced by the knowledge of the saturation scales of the colliding systems. A quantitative 
description of the multiplicity fluctuations requires more sophisticated treatment of the non-perturbative 
dynamics of multi-particle production, modeling of the details of wave functions of the collision systems, 
and fluctuations in the geometry of initial stages of collisions [siiziiHiini- 

The IP-Glasma model that combines the IP-Sat dipole model and non-linear dynamics of the Glasma 
gluon fields can provide an ab inito approach to such problems in the framework of the GGG [Bin]. 

2. DifTerent Sources of Fluctuations in CGC : the IP-Glasma model 

The input to the IP-Glasma model is the distribution of the saturation scale in the co-ordinate space 
transverse to collision direction at a given energy, inside a proton from the IP-Sat model [12]. 

The transverse distribution of Q|(s_l) is controlled by the thickness function of the proton, which in the 
IP-Sat model has a Gaussian form given by 

where the parameter Bq = 4 GeV“^ that determines the width of the gluon distribution is obtained 
from fits to HERA data mils]. The transverse distribution of Q|(s_l) inside a proton obtained this 
way is a smooth distribution. This distribution of Q 5 (s_l) can be used to sample different configurations 
of the color charges inside the proton according to the McLerran-Venugopalan (MV) model [Il[2| using a 
relation between (5|(sj_) and the variance of the color-charge distribution g^fj,{s±) in the MV model [14]. 
This is how sub-nucleonic fluctuations are introduced in the IP-Glasma model. The distributions of the 
gluon fields inside two colliding protons are obtained by solving Classical Yang-Mills equations in 2-1-1 D 
lattice. Via the solutions of these equations, the original fluctuations in the configuration of color charges 
get transmitted to the gluon helds of the two colliding protons. 

The second source of the fluctuations in the IP-Glasma model is of geometric origin. The gluon 
fields produced after the collision of two protons are obtained in terms of the gluon fields of the two 
colliding protons. The resultant field after the collision vanishes in the region where the fields of the two 
protons do not overlap. Therefore an additional source of stochasticity is introduced by the fluctuations 
in the overlap area of the two protons due to fluctuation in the impact parameter of the collision. The 
distribution of the fluctuation of the impact parameter in p-fp collisions is not known form first principles. 
However, it can be well approximated according to an eikonal model for p-fp collisions. In that approach, 
one calculates the overlap function of the two protons at a given impact parameter b = |bj^|, given by 

T^pib) = J dW r/(s^) - b^). (2) 
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Then one defines the differential probability of a p+p collision with impact parameter b as |15j 

dP 1 _ e-<^gg^gT’pp(6) 

where is the total effective partonic cross section at a given energy of collisions. The value of 

TVgCTgg is adjusted in such a way that the denominator of Eq. [^becomes the inelastic p+p collision cross 
section cr™* at that collision energy. To estimate the value of cr™* for different energies (in GeV) 
one uses the parameterization obtained from fits to the global data of total m and elastic m cross 
sections given by 

4 ^ 2 44 8 

GeV) = 25.2 + 0.05ln(-\/i) + 0.56ln^(-s/i) H-p-g H-^ mb. (4) 

^/s ' y/s ' 

This particular ansatz for implementing event-by-event fluctuations of the impact parameter in p+p 
collisions was proposed in Ref m and adopted in the IP-Glasma model m- 

The combined effect of these two sources of fluctuations leads to event-by-event fluctuations in the 
transverse distribution of the gluon fields produced by the collisions. In the IP-Glasma model, these 
gluon fields are evolved in time r according to the Classical Yang-Mills equations. The gluon multiplicity 
is calculated after an evolution time of r ^ 1/Qs when the fields becomes weakly interacting and a freely 
streaming system of gluons. The gluon multiplicity obtained in such a way fluctuates from event to event 
due to the combined effect of the sub-nucleonic and geometric fluctuations. 

The Gaussian fluctuations of the color charges in the MV model give rise to a negative binomial 
distributions (NBDs) of gluon number fluctuations for a given overlap geometry of collisions. This was 
first shown in a perturbative computation of the Glasma flux-tube picture m and was found to be 
naturally incorporated in the non-perturbative framework of the IP-Glasma model m- Due to event- 
by-event fluctuation of the impact parameter, the resultant distribution of multiplicity is a convolution 
of multiple NBDs each weighted by the corresponding values of impact parameter given by Eq.([^. 

In Ref. [iHj it was shown that these two sources of multiplicity fluctuations are not sufficient for 
explaining the width of the experimental multiplicity distribution of charged particles in p+p collisions. 
The resultant distribution from the IP-Glasma model was found to be too narrow as compared to the 
data. After introducing additional fluctuations of (5 s(s_l) according to a Gaussian distribution, the 
multiplicity distribution was found to be closer to data. As we discuss in the next section, the inclusion 
of such fluctuations has been pointed out in several previous treatments of the GGG and is the main 
topic of this paper. 
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3. Additional Sources of Fluctuations 


Incorporating event-by-event fluctuations of the saturation momentum of the proton at a given energy 
of collisions goes beyond the conventional framework of CGC. The importance of such fluctuations to 
explain the pseudo-rapidity dependence of multiplicity in p-l-A collisions has been highlighted in a very 
recent work [20]. It has been clearly demonstrated that large fluctuations in the saturation scale of 
the proton are necessary to get the right slope of the pseudo-rapidity dependence of multiplicity for 
different centralities. The intrinsic fluctuation of the saturation momentum in the CGC framework has 
also been investigated in Ref. [U] . A large increase of the intrinsic saturation scale of protons in high 

multiplicity events was shown to be essential for a quantitative description of long range two-particle ridge 
correlations in the CGC framework Emm. In was shown that in p-l-p collisions a fluctuation of proton 
saturation scale by about 5-6 times its min-bias value is required to model rare high multiplicity 
events in which the ridge like correlation was measured in the experiment. We will revisit this point in a 
later section of the manuscript. In this work we try to demonstrate the importance of such fluctuations 
for the description of the n-particle multiplicity distribution in p-fp collisions. 

The idea of intrinsic fluctuations of the saturation momentum goes back to Ref. (331 [33] , which ex¬ 
tended a correspondence between high-energy QCD and the reaction-diffusion problem that was originally 
identified in Ref. |3SJ (33] at the level of mean field approximation. The origin of such fluctuations is better 
understood within the color dipole picture In the large Nc limit the gluon distribution inside 

a target proton at a given rapidity Y can be thought of as a distribution of dipoles. This picture is 
however valid only in the dilute regime of the target, in the saturation regime the dipole picture breaks 
down. In this picture one can therefore define the scattering amplitude of the target T{r) for a dipole of 
size r. The dilute regime of the target corresponds to small values of r —>■ 0 where the distribution T(r) 
has a long tail approaching zero. In the saturation regime (large r) where the dipole picture breaks down 
T(r) approaches the unitarity limit. The saturation scale of the target at that rapidity Qs(F) is defined 
as the inverse of the dipole size r for which T{r = l/Qs{Y)) = Tg, where Ts is of the order of unity[^ 
Therefore the value of Qs is determined by the tail of the distribution T{r) which corresponds to the 
dilute regime of the target. In the dipole picture the natural variable characterizing the size of a dipole 
is expressed as ln(l/r^), therefore the tail of T(r) is actually related to the logarithm of saturation scale 

ln(g|). 

With the evolution in rapidity Y, individual dipoles inside the target start to split into new dipoles. 

^The exact value for the choice of Tg is not important, it can only introduce a logarithmic uncertainty, it is often taken 
to be 
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This nature of the splitting of individual dipoles is a purely stochastic process. The result of such splitting 
would change the distribution of the scattering amplitude for a starting value of rapidity lb given by 
T{r,Yo) to an evolved distribution at rapidity Y given by T{r,Y). 

Due to the stochastic nature of the dipole splitting, the distribution T(r, Y) will have a different shape 
if the evolution is repeated for a second time starting with the same initial distribution T(r, lb). In the 
conventional framework of CGC, the evolution of the dipole scattering amplitude is implemented by BK 
or JIMWLK renormalization equations [301 ED Ea E3 ED . However the effect of this stochastic nature 
of the dipole splitting is not incorporated in either of the two evolution equations. These two equations 
only incorporate the evolution of the average distribution of the dipole scattering amplitude (T(r))y. 
This was initially pointed out by the authors of Ref [33]. The stochastic nature of the dipole splitting 
therefore would lead to multiple distribution of T{r,Y) for a given initial distribution of the scattering 
amplitude r(r, lb)- Since the tail of such a distribution determines the effective value of the saturation 
momenta, this would lead to a distribution of ln((5|(l")) after an evolution of rapidity Y around a mean 
value of ln(((3|(y))). In reference [33] it was shown that the variance of such distribution would 
depend on steps of evolution in rapidity and the strength of coupling. The evolution of was found 
to be linear in rapidity. However the absolute value of tr^ was argued to have a non-perturbative origin. 
The effect of non-zero was conjectured to be incorporated in a Gaussian distribution describing the 
fluctuation of ln((5|) [33]. The explicit expressions of the cumulants of such fluctuations in the context of 
a stochastic reaction-diffusion model has been derived in Ref [35| . The full derivation of the distribution 
of saturation scale has been done in Ref [33] . It was shown that the leading behavior of such distribution 
is Gaussian in ln(Q|(y)) and the correction to the Gaussian fluctuations is negligible at large Y. 

In this work we introduce such event-by-event Gaussian fluctuations of the saturation momentum of 
the two colliding protons in the IP-Glasma model and study the effect on multiplicity fluctuations in p-l-p 
collisions at RHIG and LHG energies. 

4. Results 

We employ the IP-Glasma model to estimate the event by event distribution of multiplicity. The 
parameterization for the saturation scale in the IP-Sat model is obtained from the recent fits to the 
HERA data performed in Ref m- In this work we perform the simulation for p-t-p collisions at five 
different energies ^/s = 0.2, 0.9, 2.36, 7 and 13 TeV. The details of the lattice implementations of the 
IP-Glasma model can be found in Ref. [18]. The only difference from Ref. [18] is that in this work we 
use slightly finer lattices with a transverse size of L=8 fm with N=400 cells in each direction. All other 
parameters related to the lattice implementations of IP-Glasma model are kept same to what have been 
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Figure 1: Distribution of saturation scale of the proton as described by Eq. § shown for two different values of the width a. 
The tails of such distributions correspond to rare configurations of the proton with large values of saturation momentum. 


used in Ref [15] for the study of p+p collisions. 

In every event we allow the value of Q^{s±) at every point in the transverse plane to fluctuate 
according to the distribution |23| given as 


P(HQs/{Qs))) = 

yZTTa 


ln^(Q|(s_L)/(Q|(si.)}) \ 
2a2 J ■ 


(5) 


This particular choice of the distribution gives rise to a skewed distribution of Qs/{Qs) around unity. 
The distribution of Qs(sj_) for different values of cr is shown in Fig0 For simplicity, we assume a to 
be independent of the transverse distance sj_ inside the proton]^ However as described in Ref [25], we 
assume the following rapidity dependence of a 


a^{Y) = aliYo) + al{Y-Yo), 


( 6 ) 


where Y and Yq denotes the final and initial rapidity defined as V = ln(2ys/lGeV). Here we keep 
terms up to the first order in the expansion of rapidity. The contribution to a that is purely due to 
quantum evolution in rapidity is controlled by Ui, whereas the origin of <Jq is generally assumed to be 
non-perturbative. In our context both cti and erg are treated as free parameters and our goal is to 
constrain them using the experimental multiplicity distributions. 

A given distribution of Qs sets the color charge density in the IP-Glasma model. The ratio of 
saturation scale to the color charge density, is treated as a parameter which is to be 

extracted from data. A value of Qs(sj_)/g^/r(s_L) = 0.65 was previously found to provide the good fit 


^In principle a should depend on the gluon density that varies with sy 
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to the multiplicity distribution of different collision systems. However a slight discrepancy in terms of 
explaining the long tails of the probability distributions of multiplicity in case of p+p collisions was found 
in Ref [T5]. We therefore study the dependence of our results on the variation of this parameter in the 
following section and see if this parameter can be further constrained by the experimental distribution 
of multiplicity in p+p collisions. 

With this initialization we compute the gluon multiplicity per unit rapidity at rapidity y = 0 after an 
evolution time of r = 0.4 fm/c. We do not use any prescription of fragmentation and assume the gluon 
numbers to be proportional to the number of charged particles. In this work we express our results in 
terms of scaled multiplicity N/{N) to minimize the uncertainties related to the absolute normalization 
of the multiplicity. 
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Figure 2: Left: Probability distribution of scaled multiplicity in p+p collisions for 7 TeV. Experimental data points are 
from Ref. [3Zl. The two curves correspond to two different value of a parameters. Right: The same plot for p+p collisions 
at 200 GeV. The green band denotes the uncertainties due to the variation of from 0.65 to 0.45. The experimental 

data points are for the pseudo-rapidity range of |r;| < 0.5 and obtained from Ref. 13811391 . 

The minimum bias distribution of scaled multiplicity for p+p collisions at 7 TeV is shown in Figj^ 
(left). The multiplicity distribution obtained from the IP-Glasma model for two different values of a are 
shown by the histograms. It is very much evident from Figj^that in the absence of Qs-Afrctfrations the 
estimated distribution (corresponding to u = 0) is much narrower compared to the data. This result 
is consistent to the previous findings of Ref [18]. The inclusion of Qs-fluctuations with a = 0.5 leads 
to a much better description of the data up to the values of scaled multiplicity iVch/(iVch) ~ 6. For 
Nch/{Nch) > 6 we find that the using a ratio of Qs/g^y = 0.45 instead of 0.65 gives slightly better 
agreement to the data. A possible origin of such dependence on Qs/g^g can be traced back to Ref [TT] 
where the non-perturbative corrections to the Glasma flux-tube picture were found to be dominant in 
the regime of small (or large Qs/g'^y)- For large values of g^y where the Glasma flux-tube picture is 
recovered the results will be insensitive to the choice of Qs/g^g- The fact that the IP-Glasma multiplicity 
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distribution is sensitive to the choice of Qs/g'^fJ- indicates that in p+p collisions we are still in a regime 
where non-perturbative corrections are important and the experimental data need to constrain the ratio 

Qslg^g- 



Figure 3: Left: Probability distribution of scaled multiplicity in p+p collisions for 900 GeV. Right: The same plot for 
p+p collisions at 2.36 TeV. The experimental data points are for the pseudo-rapidity range of |r)| < 0.5 and obtained from 
Ref. 13711391 . 

In Figj^ (right) we show the same distribution for 200 GeV. In this case we use three different values 
of cr = 0.4, 0.45 and 0.5. We do not observe significant differences in the shape of the distributions for 
a = 0.4 and 0.45, the distribution corresponding to 0.5 slightly overestimates the data. For a = 0.4, we 
show the dependence of the distribution on the variation of the ratio Qs/g'^g- in the range of 0.45 to 0.65 
by a green band. Below Nch/(Nch) ~ 6 we see that the variation in the value of the parameter a that 
provide good fit to both the data (at 200 GeV and 7 TeV) must be less than O.I. This would restrict the 
parameter ai in Eq.([^ to be < 0.025 corresponding to the variation of rapidity Y =9.55 (^/s = 7 TeV) 
to Y =5.99 {^/s = 200 GeV). This negligible variation of af is consistent to the expectations of Ref. ^5] , 
indicating that the contribution to the a from evolution in rapidity (at small rapidity) is negligible as 
compared to the dominant non-perturbative contribution cto- A maximum variation of a from 200 GeV 
to 7 TeV to be 0.1 as suggested by Figj^ leads to a dependence like cr^(V) = 0.16 + 0.025(y — 5.99). 
This would correspond to cr ~ 0.45 at v^=900 GeV and a ~ 0.47 at y^=2.36 TeV. The corresponding 
distributions obtained for these two values of tr are shown in Fig|^ Good agreement with data again 
confirms negligible dependence of a with energy. We would however like to emphasize the fact that since 
IP-Glasma simulations are numerically intensive we do not perform fine tuning of the a parameters. 
Based on the results shown in Fig.2 and Fig.3 we would like to restrict our conclusion to the fact that a 
Gaussian fluctuation with a ~ 0.5 ± 0.1 with weak energy dependence provides a good description of the 
multiplicity distribution in p+p collisions at RHIG and LHC energies. These results indicate that for rare 














events, that are one out of ten thousand, the proton will fluctuate into a configuration which correspond 
to a saturation scale that is more than five times the magnitude of the average saturation scale. It 
is interesting to note that a similar order of fluctuation of the proton saturation scale was previously 
considered in Ref [211122] to model the high multiplicity events in which the signals of ridge correlations 
are observed in the experiments. 



Figure 4: A prediction for the probability distribution of scaled multiplicity in p+p collisions for 13 TeV. 

Finally we make a prediction for the recent 13 TeV p+p collisions at LHC using our framework. 
According to the parameterization of Eqj^ we expect the value of cr ~ 0.51 at 13 TeV. The inelastic 
cross section for 13 TeV is obtained from Eq|^ to be = 13 TeV) 76 mb. The estimated 

probability distribution for the scaled multiplicity is shown in Eigj^ In this case we have used the value 
of g^^/Qs = 0.45. The experimental data of the corrected multiplicity distribution at 13 TeV will further 
constrain the parameter cr and its energy dependence. 

5. Summary 

In this work have explored the origin of high multiplicity events in p+p collisions using the framework 
of CGC. The dominant source of multiplicity fluctuation in the framework of CGC is known to be the 
color charge fluctuations incorporated in the MV model, which with fluctuations of collision geometry 
can provide a qualitative explanation of the origin of the high multiplicity events. However, state of the 
art simulations of CGC using the IP-Glasma model demonstrate that additional sources of fluctuations 
are required to explain the origin of rare events which give rise to the tails of the multiplicity distributions 
in p+p collisions. In this work we demonstrate that including an additional source of fluctuations in the 
intrinsic saturation momenta of a proton in p+p collisions, one can very well describe the multiplicity 
distributions over a wide range of energy. In the color dipole picture of high energy scatterings in QCD, 
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such fluctuations are a consequence of stochastic splitting of dipoles that are not accounted for in the 
conventional frameworks of CGC. In this work we try to incorporate such process by fluctuating the 
intrinsic saturation momentum of two colliding protons in the IP-Glasma model according to a Gaussian 
distribution. The width of the Gaussian distribution a is an unknown parameter in this approach. We find 
that the multiplicity distributions over a wide range of energy form 200 GeV to 7 TeV is explained with a 
Gaussian of width a ~ 0.5. Two important conclusions that can be made from such observations are: (i) 
the fluctuations of 0(1) in the saturation momentum are essential to explain the origin of high multiplicity 
events in the framework of CGC and (ii) the energy dependence of such fluctuations as constrained by 
the data appears to be negligible; that is, the fluctuations are dominantly of non-perturbative origin. 
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